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This paper describes the work for optimizing the propulsive efficiency of 
flapping airfoils, i.e., improving the thrust under constraining aerodynamic 
work during the flapping flights by changing their shape and trajectory of 
motion with the unsteady discrete adjoint approach. For unsteady problems, 
it is essential to properly resolving time scales of motion under consideration 
and it must be compatible with the objective sought after. We include both 
the instantaneous and time-averaged (periodic) formulations in this study. 

For the design optimization with shape parameters or motion parameters, 
the time-averaged objective function is found to be more useful, while the 
instantaneous one is more suitable for flow control. The instantaneous 
objective function is operationally straightforward. On the other hand, the 
time-averaged objective function requires additional steps in the adjoint 
approach; the unsteady discrete adjoint equations for a periodic flow must 
be reformulated and the corresponding system of equations solved 
iteratively. We compare the design results from shape and trajectory 
optimizations and investigate the physical relevance of design variables to the 
flapping motion at on- and off-design conditions. 
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Nomenclature 


C f = Time-averaged thrust coefficient 

Cw = Power coefficient for time averaged work required to maintain the flapping motion 
c = Chord length of airfoil 

F, f= Objective functions for the sensitivity analysis and optimization 

frh = Physical plunging frequency ( =UJnc ) 

h(t)= Time dependent equation of plunging motion 

h 0 = Non-dimensional maximum amplitude of plunging motion 

k = Reduced frequency ( = TzfrhdU.-,) 

kh = Nondimensional maximum body speed ( =kho ) 

m = Index of frequency mode for superposed Fourier function 

Re = Reynolds number 

St = Strouhal number ( =2frhhoc / U x =2kho/n) 

T = Period of flapping motion 
Coo = Flight (free stream) velocity 
a = Angle of attack 

[ip = Amplitudes of superposed Fourier functions 
t] = Propulsive efficiency (=C, UJ Cw) 

T = Phase angle 

D = General form of design variable vector 
Q = Three-element vector of primitive variables (p , u, v) T . 

R = Vector of discrete residual equations for incompressible Navier-Stokes equations 
R s = Vector of discrete residual equations for Steady-State fluxes 



R t = Vector of discrete unsteady time-derivative terms 
X = Vector representing geometric variables of mesh 

I. Introduction 

D evelopment of efficient and accurate optimization approaches for aerospace applications 
is a focused interest under NASA’s current Subsonic Fixed Wing Project. A special 
challenging topic in optimization is parameterization in both spatial and temporal dimensions, 
because they directly determine the details of physics to be captured and the scales in both 
dimensions affecting the design outcome. Hence, the design of micro-air-vehicles (MAV), whose 
flapping mechanisms are inspired primarily by flying animals, not only is of considerable current 
interest technologically, but also is a unique benchmark problem for developing optimization 
methodologies. 

The interest in MAVs has become conspicuous since the 1990s [1-5]. The kinematics, 
dynamics, and flow characteristics of flapping wings have been investigated by experimental and 
computational means. Most studies have been devoted to understanding the mechanisms for lift, 
thrust and propulsive power, but only a few studies were concerned with optimization of the 
shape of flapping wing and the trajectory of motion. Lee et al. [5] designed a tadpole type of 
flapping airfoil in their parametric study based on computational fluid dynamics (CFD) analysis. 
Tuncer and Kaya [6, 7] designed a bi-plane configuration by including parameters of sinusoidal 
motion to maximize propulsive efficiency using a gradient-based approach. Most previous 
design works were conducted focusing on uniform sinusoidal plunging and pitching motions, 
where the number of design parameters is small (limited to two or three), because of a 
prohibitive computational cost in unsteady optimization problems — even with an efficient 
gradient-based method. In a realistic aerodynamic shape optimization, the number of design 



parameters needed to adequately describe an airfoil shape in transonic flow can vary from 0(10) 
to 0(100), even for a two-dimensional airfoil application [8-10]. As a result, the design cost of 
optimization based on CFD is expensive, even with today’s computation power; hence a 
judicious choice of optimization approach for the current problem is warranted. In our study, we 
adopt the approach of using adjoint sensitivity analyses [8-10] in our aerodynamic shape and 
trajectory optimizations for flapping airfoils. 

The situation in the trajectory optimization becomes more acute if one intends to explore an 
extended design space by using more non-harmonic trajectory functions. It is known that the 
flapping motion results in a strong nonlinear wake flow where propulsive thrust is generated [11, 
12]. Vytla et al. found a bifurcation behavior of vortex shedding which depends on the initial 
attitude of the airfoil [13]. In addition, they found that the bifurcation pattern can result from a 
disturbance during the flapping motion and the flapping creatures can change the flow pattern 
around their wing from one mode to a favorable one by making use of disturbances in motion. In 
other words, if the design parameters for the flapping motion include more disturbances of 
various modes, the optimal motion can be more efficient. To realize this possibility in the present 
work, the motion trajectory will be represented by considering a large number of different modes. 

The adjoint approach has been regarded as a very efficient sensitivity analysis method for 
aerodynamic shape optimization because the computational time cost is almost independent of 
the number of design variables [8-10]. To apply this method to an unsteady aerodynamic 
application, two different ways of accounting for physical time can be taken. One treats it in an 
instantaneous manner (a local approach) for dealing with phenomena in small time scale, e.g., for 
flow control. In this case, it is possible to extend the steady-state adjoint formulation in a 
straightforward fashion by adding a source-like term from the physical time difference. Kim et al. 



[8] applied this approach for the design of a high-lift device at high angles of attack, where a 
large flow separation is observed during the design optimization. The other approach treats the 
unsteady optimization problem in a time-averaged manner (a global approach), and it is used in 
this study for optimization of shape and trajectory. In order to differentiate the time-integrated 
objective function, the same number of adjoint vectors and adjoint equations as that of physical 
time steps should be prepared. Nadarajah and Jameson proposed a reverse time-integration 
method with a terminal condition to evaluate the sensitivities of the time-averaged objective 
function for two-dimensional unsteady compressible flow [14]. Their formulation, however, is 
valid only when the time integration is conducted from initial to final physical time step. Thus, 
the solution of the flow problem has to be stored for all time levels over which the optimization 
is solved. In the case of three-dimensional problems or two-dimensional problems that are very 
sensitive to time steps, the storage of flow properties can be prohibitively large. Yamaleev et al. 
[15, 16] devised a local-in-time, time-dependent adjoint method to reduce the memory size by 
dividing the whole time domain into sub-intervals for several two-dimensional unsteady design 
problems. They reformulated the sub-interval form by obtaining the terminal adjoint vector from 
previous design iteration to have a closed system equation. Their approach provides reasonable 
accuracy for the sensitivities, showing similar optimal solution to Nadarajah’s formulation. 
However, the assumption is not mathematically rigorous and may affect the accuracy of the 
sensitivity in case that flow pattern changes remarkably between the design iterations. Instead of 
assuming the terminal conditions, especially for periodic unsteady flows, we suggest to adopt 
periodic conditions of flow derivatives and adjoint vectors and to solve the sub-interval adjoint 
system equations by an iterative manner. 



In the present work, we develop an unsteady discrete adjoint code that aims to provide 
accurate sensitivities for a periodic problem such as flapping motions in a time-averaged manner. 
Based on this strategy, we perform the design for optimal shape and motion of a flapping airfoil. 
Finally efforts are made to identity parameters that are effective for increasing flapping 
performance at on- and off-design conditions. 



II. Sensitivity Analysis for Unsteady Flows 


A. Unsteady Discrete Adjoint Equations Formulated in Instantaneous Manner 

As mentioned above, an instantaneous adjoint formulation can be useful for optimal flow 
control problems by revealing the flow physics responsible for aerodynamic performances. In 
these problems, the objective function is defined as a function of flow variables at each time step, 
as in the steady-state formulation. As a result, the formulation remains essentially the same 
except that the treatment of the time derivative is included as a source-like term. The source term 
results from the physical time derivative, which in our case is approximated by a three-point 
backward Euler implicit method. As the flapping airfoil under study is moving at low speed, the 
resulting flow phenomena will be described with incompressible Navier-Stokes equations. The 
system of equations is rendered hyperbolic by incorporating the pseudo-compressibility concept 
[17]. The final discretized equations written in delta form between two consecutive iterations in a 
pseudo-time formulation with step Ax become 
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where J is the Jacobian matrix of transformation between the Cartesian and body-fitted 
coordinates, x and superscript in refer to pseudo time, and t and n to physical time. R s and R t are 
the discretized forms of the spatial-derivative and time-derivative terms in the Navier-Stokes 
equations, respectively. The discrete residual R of the unsteady flow equations is the sum of 
these two discrete derivatives, and the time-accurate solution is obtained when R = 0. For the 
design problem, the residual is a function of flow solution Q, grid position X, and design 


variables D. Hence, 



R = R S +R, =R{Q,X,D} = 0 


( 2 ) 


Similarly, the aerodynamic objective function /is also dependent on Q, X, and D, 

/ = /(Q,X,D) (3) 

Sensitivity derivatives of the aerodynamic functions with respect to design variables D are 
calculated by differentiating Eqs. (2) and (3): 

JR 3R JQ dR dX dR 

= — + + = 0 (4) 

JD dQ dD dX JD dD 

A = A^ + AA + A (5) 

dD dQ JD dX JD dD 

where the design variables vector, D, is in general time-dependent and affects the flow property 
at the corresponding physical time step. The quantities Q n _i and Q n are known and regarded as 
fixed with respect to the design parameters. Hence, we shall denote Q n +i by Q hereafter for 
simplicity. In the adjoint method, the sensitivity derivatives of the aerodynamic functions are 
obtained by combining Eqs. (4) and (5) through the introduction of the Lagrangian multiplier 
matrix A: 

df df dQ df dX df , T fdRdQ dR dX dRd 

— = — — — + — + — + A — + + — (6) 

dD dQ dD dX dD dD ^ dQ dD dX dD dD ^ 

where A is a three-element adjoint vector consisting of the Lagrangian multipliers (Xi, X 2 , Xsf 
corresponding to the primitive variables (p, u, v) T . The geometrical sensitivity vector dX/dD can 
be obtained by differentiating the grid generation code. In the present work, finite difference 
approximation is applied for simplicity. Rearranging Eq. (6) to factor out dQ/dD yields the 
following equation: 


^_^^ + _^ + Ar |"dR^ + dRY (<V_ +X T dR^ldQ 
dD dX dD dD ^ dX dD dD y ^ dQ dQ ^ dD 



By excluding the computationally intensive term dQ/dD, the sensitivity derivatives of the 
aerodynamic functions can be considerably reduced to 


# = ^ rfX + ^ +A r 
dD dX dD 3D 

This is possible if and only if the coefficient of dQ/dD vanishes for a proper vector A. After 
including the derivatives of R t with respect to Q, this becomes 

dR r 3 f dR t 1.51 

dQ dQ [ dQ JAt 

An iterative procedure will be required to obtain the solution vector A in Eq. (9), as its 
coefficient matrix can be of a large dimension. Mimicking the solution procedure for solving the 
flow variables in Eq. (1), the backward Euler implicit method with pseudo-time marching is also 
used here. Hence, 
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Note that for steady-state computations, the physical time term vanishes by setting the physical 
time step At to infinity. Boundary conditions in the adjoint variable methods are also derived 
similarly as in Eq. (9): 
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where subscript B denotes the boundary cells. 



B. Unsteady Discrete Adjoint Equations Formulated in Time- Averaged Manner 

The time-averaged formulation is useful for shape optimization where the mean performance 
is important. In this case, the objective function over a period can be expressed in a discrete 
time-integration form as 

F= £ /„(Q(0,X(0,D) (12) 

n=N m + 1 

where f is the n th physical time step within the ( m+l) th period. N m and N m +i denote the final time 
steps of the m th and (m+lf 1 periods, respectively./, is the considered function at t n (the n th time 
step) in the time averaged form for a time span of interest (7), divided by ( N m +i - N m ) time 
increments At/T. For example, for minimizing the drag coefficient, we take f n =Cd(t n ) At /T. 
Hereafter, Q(t n ) and X(t") are simply denoted by Q„ and X„, respectively. The sensitivity of the 
cost function F can be represented as 

dF = [ 9/,(Q„,X„,D) rfQ„ ( | y,(Q,,X„,D) l 

c/D n ^+i\ dQ„ dD dX n dD 3D j 


Different from the expression in the previous section, there are many flow derivative terms, 
{dQn/dD}, which are destined to be excluded in order to reduce the computation cost associated 
with these terms. Thus, these terms should be eliminated as before by combining adjoint vectors 
and residual equations. The residual equation at the n th time step with second-order time 
accuracy can be represented as 
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The sensitivity equation, after combining with the residual equations at every physical time 


step, becomes 
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If the flapping motion is completely periodic, then we can use the following identities: 
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As before, we obtain the adjoint equations for every physical time step by eliminating the 
coefficients associated with the flow derivatives term {dQ„/dD}: 
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The adjoint system equations are solved in a reverse time order by beginning with the 
equation at the ( N m +i) th physical time and ending at (N m +l) th time. Hence, this process is called 
the reverse time integration, as it traces from the later to earlier physical time steps. If the adjoint 
equations are all satisfied with the correct A„, then the sensitivity will be integrated by 
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However, the initial adjoint vectors needed to begin the time integration, A Nm +i+i and A Nm+ i+2, 
are yet to be determined. It is assumed that the initial guess on their values is immaterial and can 
be arbitrary, since they will be iteratively updated until convergence is declared. Thus, we set the 
initial values of A Nm +i+i and A Nm +i+2 to be zero and update by the periodic condition in Eq. ( 18 ) 
as the iteration goes over. The gradient can be determined if the following criterion is satisfied 
by measuring some error E between the q th and {q+\ f' steps, 

£(<,+l) = |G" , -G*|/#„< £ (20) 

where G q is the sensitivity vector at q h iteration, Ndv is the number of design variables, and £ is a 
tolerance for convergence. The boundary conditions are handled in the same manner as for the 
instantaneous formulation. 



III. Design Formulation 


A. Design Parameters 

Although a non-harmonic rigid body motion is the focus of this study, we shall use the 
uniform sinusoidal motion for plunging for validation because of availability of experimental 
data and other computational reports. The equation for a baseline rigid body motion is defined as 
follows: 


Baseline plunge motion: h = h 0 cos(2/rfr h t) (21) 

Plunge frequency: fr h = UJc / nc (22) 

where ho is the reduced plunge amplitude, k is the reduced frequency, and c is the chord length of 
airfoil. Here, several modes of sinusoidal functions with sought-after amplitudes and periods, as 
given in Eq. (23), are superimposed to the baseline plunge motion to achieve optimization. The 
general flapping motion is now represented by m additional harmonics that are compatible with 
the period of the baseline motion. 
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The total number of design parameters is 10 in the trajectory study: they are 10 amplitudes, f p . 
The constraints for fi p are limited between -1.0 andl.O, except /T whose frequency is the same as 
the baseline. Thus, we set 0 < < 2 so that the motion equation can represent non-flapping flow 

as well. Now the complete motion is expressed as follows: 
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In addition to the parameters controlling the flapping motion, the geometry of an airfoil is also 
responsible for producing thrust and lift forces, just as thin and thick airfoils will induce different 



amounts of vortex shedding. Thus, optimization for the shape of airfoil is also conducted to 
investigate its effectiveness in comparison to the flapping motion alone. The Hicks-Henne 
functions [18], as shown in Fig. 1, are adopted to describe the airfoil shape. Here 10 functions 
with different peak positions are used on the upper and lower surfaces, resulting in 20 design 
variables in total. The amplitudes of all functions, contributing to the deformation of the airfoil, 
are constrained by the limit from -0.006c to 0.006c. 

B. Cost Functions 

For the flapping motion, thrust, lift, and power coefficients are regarded as major performance 
factors. To give an overall performance through the entire unsteady motion, these coefficients 
are evaluated in a time-averaged manner. The thrust coefficient is based on the integration of the 
negative drag coefficient over a flapping period as follows: 

C^ljn-CA,^ (25) 

Here, T is the period of the motion and Q (t) represents the time-varied drag coefficient. The 
power coefficient CV is a time-averaged work required to maintain the flapping motion. 

(26) 

Here, it is assumed that the mass and inertia of the airfoil are negligible. A denotes the control 
surface of the geometry, V is the airfoil velocity, and C p is the pressure coefficient at a surface 
element with the area vector dA. The propulsive efficiency rj is defined as the ratio of the 
thrusting power to the power coefficient: 



where £/«, is the free-stream velocity or the velocity of the flapping airfoil in still air. The time- 
averaged lift coefficient is evaluated in a similar manner: 

c,=irc,(0* (28) 

1 Jt 

Depending upon the objective of an optimization problem, any one of the above performance 
indicators or a combination of them may be chosen. For a hovering motion, the lift coefficient 
and power coefficient can be considered as the objective and constraint, respectively. In the 
study presented here, the propulsive efficiency is to be maximized so that the rate of increase in 
thrust is faster than that of power input during the optimization. 



IV. Validation of Flow Analysis and Sensitivity Analysis 


A. Validation of Force Prediction 

The solver adopted for the current study is validated through predictions of time-varying 
aerodynamic coefficients on NACA 0012 airfoil for two drag inducing conditions where the 
Strouhal number is small and the Reynolds number moderate. It is reported that the contribution 
of viscous drag reaches about 70% of total drag [24] and the time-accuracy of the force 
measurement is highly sensitive to the grid distributions in these conditions. The flow is assumed 
laminar as the effects of Reynolds number did not appear significant for MAV flapping wings in 
previous work [1-7, 11-13, 20, 21]. We specify the flapping motion by imposing the velocity 
and displacement at grid points on the airfoil according to Eq. (24). Since there is no 
experimental data provided, we compared the current force measurement with the previously 
reported CFD results by Young [23] and Lian et al. [24] The first case is at the reduced 
frequency, k, of 3.93, non-dimensional plunging amplitude, ho, of 0.0125 so that the resulted 
Strouhal number, St, is 0.03. Figure 2 presents the results of Q and Q from three different 
densities of C-type grid: coarse (257x49), intermediate (513x49), and fine (1025x193) and a fine 
(513x173) O-type grid. Besides the coarse C-type grid where the number of grid points along the 
airfoil direction is not enough to provide an accurate prediction of drag force, solutions of all 
other grids fall between those reported in [24], The lift coefficients of all grids, including coarse 
grid, show high coincidence with each other. Figure 3 shows the histories of drag and lift forces 
at the same reduced frequency but with a doubled plunging amplitude (hence doubling St) of 
ho= 0.025, and Strouhal number, St, of 0.06. In this condition, the discrepancy between coarse 
grid and other grids is considerably reduced from the previous case. Thus, it can be reasoned that 
the effect of viscous force gets less dominant and the grid dependency gets weaker as the 



Strouhal number increases. Since the purpose of current study is to perform design optimization 
for propulsive efficiency, the force prediction is validated over a wide range of conditions 
covering the most-studied flapping flight regimes. The efficiency curve with respect to kh (non- 
dimensional flapping velocity) in Figure 4 confirms the force prediction of current study is well 
validated. Here, the non-dimensional body speed is usually calculated by a multiplication of 
reduced frequency and wake width, which is the excursion distance of the trailing edge. For the 
pure plunging motion, the plunging amplitude is regarded as the excursion distance. The reduced 
frequency is usually calculated by InfruclU <». However, in the present work it is defined as 
nfrhdUn, in order to be consistent with that used in the experimental and other CFD results [13, 
20, 21]. Thus, the kh value in the present work is half of that in Ref. [23]. We compare the force 
predictions on three C-type grids with different sizes, 513x49, 513x97 and 1025x193, with 
experimental data [11] and Young’s CFD computation [23] at the range of kh from 0 to 0.6 
(0-1.2 in [11, 23]) at the Reynolds number of 20,000. The plunging amplitude is fixed at 0.175 
and reduced frequency is varying from 0 to 3.4. Similarly to results in previous reports, the peak 
of efficiency is obtained around k/?=0.1875 (0.375 in [11, 23]). In the following Section V, an 
optimization case of symmetric airfoil will be presented at this highest efficiency point. 

B. Validation of Flow Analysis for Wake Profiles 

In order to study the design of a more practical case of lifting airfoil, the asymmetric SD7003 
airfoil [19] was used in our simulations, while previous experimental works [13, 20, 21] were 
used as basis for comparison. The Reynolds number for all validation cases of plunging motion 
is 10,000. The reduced frequency, k, is 3.93; the plunging amplitude, ho, is 0.05; and the angle of 
attack is 4°. An O-type structured grid around the airfoil is used with three grid sizes of 243x126 
(coarse), 243x265 (medium), and 486x358 (fine) for the convergence study. For validation of 



the baseline flow solver, the CFD results are compared with experimental data for a pure plunge 
motion [19]. 

The computed vorticity contours at selected physical time steps for the plunge motion is 
presented in Fig. 5(a), in which the CFD results from the 20 th cycle of the motion on the fine grid 
are used. The experimental contours, shown in Fig. 5(b), are the phase-averaged particle image 
velocimetry (PIV) data over 120 periods. It shows that the computed solutions at various phases 
are in good qualitative agreement with the experimental observation. The wake profiles at two 
chords away from the leading edge are compared in Figs. 6, for each of the four phases given in 
Figs. 6. In addition, the results of three grids are included, indicating that the flow patterns are 
essentially grid converged, except at V F = 270°, where the fine-grid solution departs from the 
other two solutions significantly; overall, even the coarse-grid solution appears to capture 
essential physics well. Hence, in the interest of reducing computation time, we chose to use the 
coarse grid for the optimization. Furthermore, the resolution of the physical time step is another 
factor that can affect the accuracy of an unsteady flow analysis. The time steps selected are 60, 
100, 200, 500, and 1000 steps for a period, and the corresponding results are given in Fig. 7. The 
profile predicted by 200 physical time steps is nearly identical to that of 1000 physical time steps. 
However, the results of 100 and 60 steps show the structure of wake vortices as somewhat 
dissipated. Hence, 200 physical time steps in a period will be used in the design optimization. 

C. Validation of Sensitivity Analysis for Steady and Unsteady Incompressible Viscous 
Flows 

Next, we present the validation of the sensitivity analysis, based on the discrete adjoint 
approach, for steady as well as unsteady flows. For both cases, geometrical variables are taken as 
design variables. In order to validate the differentiation of the incompressible RANS solver, the 



sensitivities from the adjoint solver are compared with finite-differenced gradients (referred to as 
FDM) for a steady flow over the SD7003 airfoil at the angle of attack of 5° and Reynolds 
number of 10,000. The objective function is the lift coefficient, and 20 Hicks-Henne functions 
are used as design variables (10 on the upper and 10 on the lower surfaces, respectively). As 
shown in Fig. 8, all the sensitivities by the adjoint solver are in very good agreement with the 
results by FDM with a step size of 0(1 O' 7 ). 

For the sensitivity analysis of unsteady flow, the same geometry is applied at an angle of 
attack of 4° as the initial posture, with a plunging amplitude of 0.05, reduced frequency of 3.93, 
and Reynolds number of 10,000. As mentioned above, in order to prevent possible violation of 
periodic assumption of current adjoint formulation during the flow and sensitivity analyses, only 
the plunge motion, the simplest flapping motion, is considered for the validation and design 
optimizations of the present work. It is reported that a nonlinear bifurcation in the wake flow can 
happen during some pitching or pitching combined with plunging motions at certain conditions 
[13]. The sensitivities from the time-averaged approach are compared with the FDM with the 
step size of 0(1 0' 5 ) in Figs. 9-11. Recall that the motion and flow pattern are assumed periodic 
in the present adjoint formulation. Therefore, the time averaged coefficients at the 20 th flapping 
cycle, at which the motion has attained periodic characteristics, is used for the validation and 
design optimization. Figure 9 is the sensitivities of lift coefficient with respect to the shape 
parameters. The design variables numbered from 1 to 10 denote the geometric changes from 
leading edge to trailing edge on the upper surface, and those of 1 1 to 20 refer to the geometric 
changes on the lower surface. It is shown that the sensitivities show a good comparison with the 
FDM (step size of 0(1 0' 5 )) except in the upper leading edge regions where complicated vortex 
structures are observed. In other words, the region is still showing aperiodic characteristics, and 



the assumption of periodic condition for the adjoint formulation is not satisfied locally. To 
investigate the effect of local aperiodic flow characteristics on the sensitivities, we present the 
history of thrust coefficient and its finite differenced sensitivities with respect to select design 
parameters during 20 th cycles in figure 10. It can be observed that the FDM sensitivities are still 
varying considerably even between 19 th and 20 th cycles where the thrust coefficient already 
converges to a certain value with a very small, i.e., O(10' 5 ), amount of deviation. Among the 
design parameters, the leading edge parameters on the upper surface denoted by parameters 1, 2, 
4 show larger variations while others, 5 and 8, which are relatively closer to the trailing edge 
region, have already converged. Thus, it can be reasoned that a small local aperiodic flow pattern 
(possibly due to both of the numerical or physical aperiodicity) affects the accuracy of the 
sensitivities, especially at the leading edge region. Overall, however, the direction of the gradient 
vector is almost same as those calculated by FDM. The convergence characteristics of the 
sensitivities of using different numbers of outer iterations are presented in Fig. 11. In the figure, 
the adjoint sensitivities of the thrust coefficient with respect to the shape parameters from each of 
outer iterations are compared with FDM. The first iteration with the initial conditions A Wm +i=0 
and A, V/ „ + 2=0 clearly provides a poor accuracy. As the number of iterations increases, the 
sensitivities converge rapidly. After the 4 th iteration, the accuracy of sensitivities appears nearly 
unchanged if E < 0.01. Hence, the criterion E < 0.0 1 is used throughout the study. 



V. Optimization and Performance at On- and Off-design Conditions 


For the design of a flapping wing with plunging motions, two shape optimization cases and 
one trajectory optimization are considered to demonstrate the capability of the current design 
approach. The objective of all optimization cases is to maximize propulsive efficiency. Shape 
optimizations are carried out on symmetric and asymmetric airfoils, respectively, by using 20 
Hicks-Henne functions. The constraints include bounds on the geometrical parameters to prevent 
too severe change in shape. For the trajectory optimization, non-harmonic motion is sought after 
by using the amplitudes of 10 superposed Fourier functions as design parameters. The frequency 
and amplitude of mean motion are not included in the set of design variables because those 
effects will be investigated for off-design conditions by varying kh, a representative non- 
dimensional body speed. The optimization algorithm used is the BFGS (Broyden-Fletcher- 
Goldfarb-Shanno) non-constrained line search method [22], 

A. Shape Optimization for a Symmetric Airfoil (NACA 0012) 

As a validation of shape optimization capability, a shape design based on NACA 0012 is 
conducted to obtain the optimal shape in view of propulsive efficiency. The design conditions 
are: Re = 20,000, k = 1.07125, ho = 0.175, and a = 0°. As observed in Fig. 4, these conditions 
correspond to the highest peak of propulsive efficiency of the baseline NACA 0012, the shape 
optimization will attempt to further improve the maximum efficiency. The intermediate C-type 
grid (513x49) is used for the optimization, in view of balance between the computational 
accuracy and cost. Then, the performance and contours are obtained by using fine grid to 
evaluate improvement of the optimized design. Figure 12 shows the design history. The 
efficiency converged to 0.3157, an increase by 11% from 0.2844, while the thrust increases by 
8.3% from 0.06205 to 0.06767 and power coefficient reduces by 1.8% from 0.2182 to 0.2143, 



after about 33 design iterations. The profile of the optimized airfoil is compared with the NACA 
0012 in Fig. 13. The key features of the geometrical change in the optimized airfoil are: the 
leading edge grows thicker and the trailing edge thinner than the baseline airfoil. It is noticeable 
that the resulted geometry is very similar to a tadpole-like airfoil, as suggested by Lee et al. [5]. 
In their parametric study, it was found that a thick leading edge is favorable for better propulsive 
efficiency in a relatively high-frequency flapping motion because a thick leading edge can 
reduce the strength of vortices generated there, which increases required power input. Regarding 
the trailing-edge shape, they suggested that a thinner airfoil could produce a larger pressure 
gradient between upper and lower surface, which in turn can generate stronger thrusting jets 
responsible for producing a positive net thrust force. The time histories of thrust coefficients by 
the baseline and optimized airfoils are compared in Fig. 14. Figure 14 shows that the symmetric 
geometry results in an almost identical thrust generation at up and down strokes (i.e., symmetry 
in C t ). The improvement performance by the optimized geometry is eminent at positive peaks in 
the thrust coefficient (e.g., about t=9.25T and 9.75T). Figures 15-17 suggest how the wake of the 
optimized airfoil differs from that of the NACA 0012. In the current design condition, the 
leading edge vortices (LEVs) flow down and interact with the trailing edge vortices (TEYs). It is 
observed that the LEV, as it convects downstream towards the trailing, leads to the formation of 
vortex pairs. For both airfoils, the wake vortex structures, after detaching from the trailing edge, 
show complicated patterns attributable to the LEVs. However, they eventually evolve to form the 
inverse Karman vortices as the small vortex singlet and/or pairs are merged or dissipated as 
shown in Fig. 15. For a clear definition of wake trend, we shall use “S” to denote a single vortex 
and “P” to denote a pair of vortices of opposite signs [25]. Accordingly, in Fig. 16, the baseline 
airfoil shows a 6P+2S wake (3P+S at each stroke) pattern at the moment right after the vortices 



are detached from the trailing edge. On the other hand, the optimized airfoil shows a 10P+2S 
wake pattern (5P+S at each stroke) temporarily. Consequently, the generation of 2 extra vortex 
pairs per stroke by the optimized airfoil enhances its thrusting performance. 

The mechanism of improvement can be explained by the role of LEYs. On the baseline airfoil, 
two LEVs (A and B) flow down to the trailing edge at t/T=9.5 on the upper surface and at 
t/T=10.0 on the lower surface (C and D), respectively. The LEVs, A and C, originate in the very 
leading edge region and the others, B and D, near x/c=0.25, where the airfoil has maximum 
thickness, at around t/T=9.35 on the upper surface (see Fig. 15) and at about t/T=9.85 on the 
lower surface. With the shape change, the optimized airfoil gets thicker in the front section and 
the trailing edge becomes sharper. As observed in Figs. 15 and 16, the LEVs in B (see t/T=9.5) 
and D (t/T=10) become larger than those of baseline airfoil, while the LEVs in A (t/T=9.75) and 
C (t/T=9.25) get smaller. The strengthened LEVs B and D make a ‘2P’ in the wake by 
interacting with TEVs, instead of a stretched ‘P’ in the baseline airfoil. As shown in Figs. 17, the 
vorticity contours and streamlines at t/T=9.3 indicate that ‘2P’ structure of optimized airfoil is 
much more favorable than ‘P’ wake pattern of NACA 0012; the arrows in the figure indicate the 
velocity collectively generated by the pair at the neighboring boundary. It shows that addition of 
positive velocity has been achieved by the optimized airfoil, rather than cancellation in the case 
of NACA 0012. Thus, the effect of thrusting jet is maximized on the optimized airfoil by 
synchronizing generation of vortex pairs by the shape with the imposed flapping motion to yield 
a favorable wake field. Here, we define the favorable wake field as having thrust-inducing jet 
flows in which a counter-clockwise vortex is placed above the clockwise vortex so that the 
induced jet flow is directed towards the stream direction. 



B. Shape Optimization for an Asymmetric Airfoil (SD7003) 

Now we shall turn our attention to the designs of an asymmetric airfoil, SD7003, for a lift 
generating flight. The design conditions for shape and trajectory optimizations remain the same 
as those for validation study in Section IV-C: Re = 10,000, k = 3.93, ho = 0.05, but a = 8° 
allowing a higher lift force. The design history for maximizing propulsive efficiency with shape 
parameters is presented in Fig. 18. In less than 20 design iterations, the process has converged; 
the propulsive efficiency is increased by 14.2% from 0.1584 to 0.1809, and the thrust coefficient 
improved by 13.4% from 0.04605 to 0.05221. Both performance indicators have very similar 
convergence characteristics because the power input has kept within a very slight variation from 
the baseline geometry, decreasing from 0.2907 to 0.2886. For a further investigation of geometry, 
the airfoil profiles of the baseline and optimized geometries are compared in Fig. 19. The overall 
trend for shape change is very similar to those of previous case, i.e., a thickened leading edge 
and a sharpened trailing edge. In addition, the optimized airfoil shows a larger change on the 
lower surface than on the upper surface, resulting in a tadpole like shape. Figure 20 compares 
the time-histories of thrust coefficients of the optimized airfoil with that of SD7003 between 
t/T=18.0 and 20.0 (19 th and 20 th cycles). Besides the interval between t/T=l 9.4-1 9.75 
(At/T=35%), the overall thrust generation increases. Especially, the thrust during down-stroke 
(t/T=19~19.4) is improved, primarily resulting from the change on the lower surface. 

The unsteady patterns of vortices around the optimized airfoil at select physical time steps are 
compared with those of the baseline airfoil in Fig. 21. Overall, both of the airfoils produce 
similar inverse Karman vortices wake pattern but the spaces between vortex cores (in both 
horizontal and vertical directions), which are generated at the same time interval, are wider on 
the optimized airfoil. Thus, it can be confirmed that the wake flow acceleration of the optimized 



airfoil is stronger than that of the baseline. A close examination of Fig. 21 reveals that the 
strength of the LEV on the upper surface of the optimized airfoil is smaller than that of the 
baseline model. This reduces power input needed for generating those vortices. Besides, the TEV 
of the optimized airfoil, which are detached from the airfoil at the end of strokes (t/T=19.5, 
t/T=20.0), are more stretched with stronger intensity. As a result, the flow acceleration could be 
enhanced by the shape optimization. The x-momentum contours plotted at t/T=19.3 in Fig. 22 
show that the optimized airfoil produces larger areas of high velocity in the wake vortices than 
the baseline airfoil. Further investigation about the optimized airfoil and the effect of trajectory 
design will be discussed in the following section. 

C. Trajectory Optimization (SD7003) 

Trajectory and geometric parameters can be considered together in the design of an MAV, 
where an instantaneous optimal combination of both parameters during the flight will be sought 
after. However, it is instructive to carry out separate optimizations for geometry and trajectory so 
that insight into individual effects can be gained, prior to taking both into account. In our 
trajectory optimization, the reduced frequency and plunging amplitudes will be fixed for each 
optimization; but they will be varied when off-design performances are studied. 

Through the trajectory optimization, the thrust coefficient increases from 0.04605 to 0.08099, 
a remarkable 75.8% increase with respect to the baseline, as displayed in Fig. 23. However, the 
propulsive efficiency increases by 7.3% only, from 0.1586 to 0.1701, because the change of 
motion parameters increases not only the thrust force but also the power input. The final design 
parameters after optimization are presented in Table 1. The design parameters that have larger 
changes than the others are selected for presentation. The major change comes from not only low 
frequency (m= 1), which is same as the baseline motion, but also from higher frequencies (m = 2 



and 5), which lead to a more oscillatory motion. Figure 24 compares trajectories of the baseline 
and optimized motions; the optimized trajectory is noted for its asymmetrical behavior, while 
modulating the baseline motion with small amplitude oscillations of higher harmonics. The 
effect of higher harmonics is more noticeable during down stroke. The vorticity contours at 
select physical time steps are compared in Fig. 25. Even though the wavy motion generates many 
small vortex pairs (‘P’) or singlets (‘S’) during a cycle, it still results in ‘2S’ inverse Karman 
vortices similarly to the baseline harmonic motion overall. However, the distance between 
positive and negative vortex cores is wider and the size of vortices is much larger than those in 
the harmonic motion. The oscillatory motion spread the small vortices in a wider wake area, 
which eventually merge together to form a larger ‘2S’ wake. Investigating closely, several jet 
producing vortex pairs are generated even in one excursion from top to bottom of the motion, 
especially at the phase from t/T=19.25 to 19.75. Those vortices are shed at the phase from 
t/T=19.75 to 20.25 (same as 19.25). Thus, most of the thrust force is generated at the phase from 
t/T = 19.25 to 19.75, so is in the baseline motion, but the wavy motion generates the thrust more 
intermittently, thereby producing more of thrust-favoring vortices. This mechanism is more 
clearly observed in Fig. 26, in which the time-varying and time-averaged thrust coefficients at 
sub-intervals for both the baseline and optimized motions are shown. The sub-interval of time- 
averaged values is taken to be 0.1T, corresponding to the smallest period (highest frequency 
mode). The added high-frequency components increase the overall amplitude of thrust 
coefficient. The maximum thrust coefficient reaches about 3.4 instantly, about 10 times that 
(0.35) of the baseline motion. While the wavy motion causes large-amplitude fluctuations in 
thrust, the end result is that a significant thrust generation over a period. As clearly revealed in 
Fig. 26, the major gain in thrust happens during t/T = 19.10-19.40, especially 19.30-19.40, in 



which the peak value in the optimized motion is more than doubled. It is noted that the thrust 
generation occurs during the down stroke phase (t/T=19.0-19.5) and during the beginning of up 
stroke (19.5-19.75) in both baseline and optimized motions. The thrust generating mechanism of 
the optimized motion can be explored by investigating the vortex structures around the trailing 
edge at select peak of thrust during this positive net thrust phase, as displayed in Figs. 27. The 
baseline motion generates ‘2S’ wakes consistently while those of the optimized motion evolve 
dynamically. Also, the streamline patterns in the optimized motion, reflecting the vortex strength 
and direction, suggests a wider jet flow area than the baseline motion does. Furthermore, small 
‘P’s and ‘S’s in the trailing edge region enhances the wake velocity components in the 
streamwise direction, especially at t/T=19.125 and 19.325, when two highest peaks of thrust are 
observed. To give a succinct comparison of effectiveness of applying optimization, the time- 
averaged performances of the shape-optimized, trajectory-optimized, and baseline models are 
presented in Table 2. The trajectory optimization doubles the thrust of the baseline motion and is 
significantly higher than the shape optimization; its lift force however deteriorates by about 1.7% 
and 2.5% in comparison with the baseline and shape optimization respectively. To sustain the 
motion, the power input needs to be increased, thereby resulting in a not as high an increase in 
thrust, but still a significant improvement, 11% and 7.3% respectively by the shape and 
trajectory optimizations. 

As impressive as it may be, it begs the question as to how well the optimized configurations 
(airfoil shape and trajectory) described above will perform at off-design conditions, as shall be 


discussed in what follows. 



D. Off-design conditions 

The performances of optimized geometry and trajectory are investigated at various values of 
kh by varying reduced frequency k from 1.965 to 15.75 (k= 3.93 is the design condition). Figures 
28 show the overall efficiency and thrust performance vs. kh for four cases: a combination of 
baseline and optimized shape and trajectory; all are computed at the same angle of attack of 8° 
and Reynolds number of 10,000. Recall that the baseline airfoil is denoted as SD7003. It is noted 
that the in the fourth case, it is assumed that the optimal trajectory obtained for the baseline 
airfoil maintains its high performance for the optimized airfoil. 

In Fig. 28-(a), the optimized airfoil gives a superior performance to the baseline airfoil at 
almost all of the conditions. Especially, the maximum efficiency is improved from 18.45% at kh 
= 0.2260 (baseline) to 22.1% at kh = 0.2555(optimized); the optimized airfoil surpasses the 
baseline in efficiency in the entire range of kh. The thrust coefficients of the optimized geometry 
are also higher than those of baseline except when 0.4 < kh < 0.55 where the SD7003 airfoil 
shows unusual pattern in efficiency and thrust coefficient. It is reported that the efficiency curves 
of flapping airfoils have only one extrema in the kh space [23], but SD7003 has two; this 
behavior is correlated with that in the thrust coefficient curve, in which it exhibits a reflection 
point. This peculiarity disappears in the optimized airfoil. Moreover, the optimized airfoil in 
baseline motion shows a superior performance to that by the baseline airfoil in optimized 
trajectory, except in the low flapping speed regime, 0.14 <kh< 0.21, as evident in the amplified 
view displayed in Fig. 28-(b). In this figure, the best performer belongs to the optimized shape 
with the optimized trajectory. These results suggest that the trajectory optimization should be 
carefully carried out via a more rigorous design formulation in which both the shape and 
trajectory optimizations are carried out simultaneously. That is, the interaction between 



geometrical and dynamical (via trajectory of motion) factors may not be negligible if a more 
robust optimization is sought after. Nevertheless, the present study already shows that the 
geometric optimization obtained at one design point is shown to be sufficiently robust that its 
superiority in performance persists in all off-design points. 

VI. Conclusion 

Optimal designs for airfoil shape and flight trajectory under periodic flapping motion have 
been conducted by using the unsteady discrete adjoint approach in which time-averaged 
objective functions are considered with periodic conditions applied to flow derivatives and 
adjoint vectors. Accurate solutions of sensitivities have been obtained with this approach in an 
iterative manner. With the adjoint approach described here, optimization of unsteady problems 
that contain a large number of design variables can be solved within a practically acceptable time 
frame. For flapping airfoil problems, improvements for the thrust and propulsive efficiency were 
achieved by shaping airfoil and/or changing the motion trajectory. 

Design optimizations are conducted for shape designs of symmetric (NACA 0012) and 
asymmetric (SD7003) airfoils and a trajectory design. NACA 0012 is re-designed at its best 
efficiency point (kh=0. \ 87) to improve the maximum efficiency. The propulsive efficiency is 
improved by 11% and the thrust increases by 8.3% and the power coefficient reduced by 1.8%. 
The improvement in the propulsive efficiency comes from higher thrust and lower aerodynamic 
work; the former originates from the sharpened trailing edge and the latter from the thickened 
leading edge.. 

For further practical applications, i.e., optimization of lift-generating airfoil shape for 
plunging motion and its trajectory are carried out based on SD7003. Trajectory optimization 
produces an increase in thrust by 76%. As more power input to sustain the maneuvering of the 



optimized motion is required, the improvement in the propulsive efficiency is about 7.8%. On 
the other hand, although the shape optimization increases the thrust force by just 13.4%, it 
produces a higher propulsive efficiency than the trajectory optimization because the power input 
remains almost unchanged with the shape. In addition, the study of off-design conditions showed 
the shape-optimized airfoil is quite robust that its performance is almost always superior to that 
of the baseline airfoil. This finding suggests that shape optimization could be very useful for 
designing aerodynamically efficient bio-inspired MAVs, when constrained by conditions, such 
as limit of power input or trajectory prescribed by an actuator or given mission. In the optimal 
design of trajectory, it is observed that the effect of oscillatory non-harmonic motion is to shift 
the frequency of motion instantaneously so as to find a locally better performance. 
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Tables 


Table 1 Motion parameters for trajectory-optimized case 

(Select sinusoidal parameters, design condition: Re = 10,000, k = 3.93, h 0 = 0.05, a = 8°) 


Parameters 

Index of design 

variable 

Baseline 

Trajectory 

optimized 

/?isin(27i/iv) 

1 

0.0 

5.1736E-03 

/? 2 cos (2nfr h t) 

2 

0.0 

-8.3726E-04 

fj3sm(4nfr h t) 

3 

0.0 

-1.4634E-03 

fpcos(4n fri,t) 

4 

0.0 

1.1642E-03 

/? 9 sin( 1 0nfr/ t t) 

9 

0.0 

-1.4559E-03 

/3iocos(l0nfri,t) 

10 

0.0 

1.9877E-03 


Table 2 Comparison of performances 

(Design condition: Re = 10,000, k = 3.93, h 0 = 0.05, a = 8°) 


Performance 

Baseline 

Shape optimized 

Trajectory 

optimized 

Q 

4.6047E-02 

5.2213E-02 

8.0991E-02 

c, 

7.0569E-01 

7.1178E-01 

6.9392E-01 

rj 

15.84% 

18.09% 

17.01% 
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Fig. 1. Hicks-Henne functions for shape optimization 
[same 10 functions are imposed on upper and lower surfaces] 
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Fig. 2. Drag and lift coefficient history 
[Re = 20,000, k = 3.93, h 0 = 0.0125, a = 0°, NACA 0012, harmonic motion] 
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Fig. 3. Drag and lift coefficient history 


[Re = 20,000, k = 3.93, h 0 = 0.025, a = 0°, NACA 0012, harmonic motion] 
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Fig. 4. Flow analysis and validation for pure plunging motion. Comparison of propulsive efficiency at 
various flapping conditions with experimental and numerical results for the time averaged performance. 
[NACA 0012, C-type, Re = 20,000, h 0 = 0.175, k= nfr h clU m , kh corresponds to half of those in [1 1,23]] 
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Fig. 5. Comparison of vorticity contours between (a) CFD and (b) experiment for pure plunging motion 
[Re = 10,000, k = 3.93, h 0 = 0.05, a = 4°, harmonic motion, CFD: 20 th cycle, fine grid, experiment: phase-averaged 

PIV for 120 periods] 
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Fig. 6. Grid refinement study of wake profile validation 
[Re = 10,000, k = 3.93, /z 0 = 0.05, a = 4°, harmonic motion, 
experiment: phase-averaged PIV for 120 periods, CFD: 20 th cycle, fine grid ] 
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Fig. 7. Parametric study of wake profile validation based on the physical time step size 
[Re = 10,000, k = 3.93, h 0 = 0.05, a = 4°, harmonic motion, 
experiment: phase-averaged PIV for 120 periods, CFD: 20 th cycle, fine grid ] 
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Fig. 8. Sensitivity analysis for lift coefficient with respect to shape parameters. 
[SD7003, ite=10,000, a = 5°] 
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Fig. 9. Sensitivity analysis for time-averaged lift coefficient with respect to shape parameters (s = 10' 4 ). 
[SD7003, ite=10,000, k =3.93, h 0 = 0.05, a = 4°] 



Fig. 10. Convergence history of time-averaged thrust coefficient and its finite differenced sensitivities 
with respect to shape parameters. [SD7003, Re=\ 0,000, k =3.93, h 0 = 0.05, a = 4°] 



Fig. 11. Convergence of the sensitivities according to the number of outer iteration 
(time-averaged thrust coefficient, shape parameters). 

[SD7003, ite=10,000, k =3.93, h 0 = 0.05, a = 4°] 



Fig. 12. Design history for maximizing propulsive efficiency with shape parameters. 
[Re = 20,000, k = 1.0712, h 0 = 0.175, a = 0°, Baseline: NACA 0012] 




Fig. 13. Design history for maximizing propulsive efficiency with shape parameters. 
[Re = 20,000, k = 1.0712, h 0 = 0.175, cc = 0°] 



Fig. 14. The comparison of time-history of thrust coefficient [two periods, Re = 
20,000, k = 1.0712, /z 0 = 0.175, a = 0°, shape optimization, fine grid] 



(b) Optimized 

Fig. 15. Comparison of vorticity contours between baseline and optimized geometry for harmonic motion 
[Re = 20,000, k = 1.0712, h 0 = 0.175, a = 0°, 10 th cycle (t/T=9.325), shape optimization, fine grid] 
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Fig. 16. Comparison of vorticity contours between baseline and optimized geometry for harmonic motion 
[Re = 20,000, k = 1.0712, h 0 = 0.175, a = 0°, 10 th cycle, shape optimization, fine grid] 











Fig. 17. Comparison of vorticity contours and streamlines between baseline and optimized geometry for harmonic 
motion [Re = 20,000, k = 1.0712, ho = 0.175, a = 0°, t/T=9.3, shape optimization, fine grid] 


0.185 


Thrust Coefficient 
Objective Function 


-,0.06 


0.155 



- 0.055 


- 0.05 O 


10 15 

Design Iterations 


20 


0.04 


Fig. 18. Design history for maximizing propulsive efficiency with shape parameters. 
[Re = 10,000, k = 3.93, h 0 = 0.05, a = 8°, Baseline: SD7003] 




(b) Shape of airfoils 

Fig. 19. Comparison of airfoil shape between baseline and optimized geometries. 
[Re = 10,000, k = 3.93, h 0 = 0.05, a = 8°, Baseline: SD7003] 
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Fig. 20. Comparison of vorticity contours between baseline and optimized geometry for harmonic motion 
[Re = 10,000, k = 3.93, h 0 = 0.05, a = 8°, 19 th and 20 th cycles, shape optimization, fine grid] 
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Fig. 21. Comparison of vorticity contours between baseline and optimized geometry for harmonic motion 
[Re = 10,000, k = 3.93, h 0 = 0.05, a = 8°, 20 th cycle, shape optimization, fine grid] 
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Fig. 22. Comparison of X-momentum contour between baseline and optimized geometry for harmonic motion 
[Re = 10,000, k = 3.93, h 0 = 0.05, a = 8°, t/T=19.3 cycle, shape optimization, fine grid] 
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Fig. 23. Design history for maximizing propulsive efficiency 
with motion parameters. [Re = 10,000, k = 3.93, h 0 = 0.05, a = 8°] 



t/T 


Fig. 24. Comparison of trajectory for two periods of flapping motion. 
[Re = 10,000, k = 3.93, h 0 = 0.05, a = 8°] 
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Fig. 25. Comparison of vorticity contours between (a) baseline and (b) optimized geometry 
[Re = 10,000, k = 3.93, h 0 = 0.05, a = 8°, 20 th cycle, trajectory optimization, fine grid] 
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Fig. 26. Comparison of aerodynamic coefficients according to physical time step. 

The sub-intervals for the time-averaged values are 0.17, which is minimum period 
among the superposed sinusoidal functions [Trajectory optimization, Re = 10,000, k = 3.93, a = 8°] 
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Fig. 27. Comparison of vorticity contours and streamlines between baseline and optimized trajectories at trailing 

edge during the positive net thrust generation 
[Re = 10,000, k = 3.93, h 0 = 0.05, a = 8°, 20 th cycle, trajectory optimization, fine grid] 
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(a) Off-design performance 



(b) Off-design performance around the design point 
Fig. 28. Comparison of aerodynamic performance at on/off-design conditions, (b) showing the expanded 
view near the design point [Re = 10,000, k = 3.93, a = 8°] 


